maybe rather put in src/sage/combinat/q_analogues.py ?

Pardon my ignorance: Why?

because similar polynomials in q (cardinal of finite field) are already there..

and your method should take q as a parameter, and work with q an indeterminate

I have a few of comments:

1. why not use just one-liner:

return sum( moebius(n//d) * q**d for d in n.divisors() ) // n

instead of the current:

    r = ZZ.zero()
    for d in n.divisors():
        r += moebius(n//d) * q**d
    return r // n

2. The documentation should state that polynomials are monic.

3. There exists a formula for the number of monic irreducible multivariate polynomials. It can be seen here or in this paper. I'd suggest to implement it as well, and add the number of variables as yet another function argument with the default value 1.
   You may find my implementation of this formula in PARI/GP at this link

Thanks for the suggestions: All done, except for the multivariate case, which is left for future work.

Depending on how important the multivariate case is, this could be modified to include m = 1 as an optional parameter with a "NotImplementedError" with comments pointing to the references Max included. Personally I think this can be left for someone else to do though?

One thing I didn't understand was whether the formula in Bodin's work (which is only asymptotic) was exactly the same as Max's comment in the forum, or whether Max's comment and PARI implementation were precise?

Personally I think this univariate one is nice as it is, especially with the polynomial output when q=None and this recursive (potentially only asymptotic) formula for the multivariate one can be solved by someone else when they're interested in doing so?

One thing I didn't understand was whether the formula in Bodin's work (which is only asymptotic) was exactly the same as Max's comment in the forum, or whether Max's comment and PARI implementation were precise?

@GiacomoPope: Bodin's result is not only asymptotic, Lemma 4 in his work gives an exact formula. That formula looks similar to mine (which was posted in the dxdy.ru forum). Historically, my formula appeared first (in 2006), but I considered it minor and did not bother to publish it elsewhere. PARI/GP code implements my formula and it was used to produce terms for many related sequences in the OEIS, some of which (such as A115457) I added back in 2006. I hope this explains the context.

Btw, A115457 links another relevant paper Survey on counting special types of polynomials.

Thanks for the context Max. This helps. Potentially this function would be useful for people, but I don't think it needs to be included into this PR for this contribution to have value. It can be added at a later date without adding any additional work.

If you disagree though, then maybe porting your PARI code for this PR is something that needs to be done.

I'm not familiar enough with Sage development, and thus I'd appreciate it if someone can properly incorporate my code into this PR. Please find my extension to arbitrary positive m below. In the original code I've changed:

• q = ZZ['q'].gen() to q = QQ['q'].gen() since in general for m > 1 we may have an integer-valued polynomial but with rational coefficients
• ZZ(n).divisors() to n.divisors() because conversion n = ZZ(n) is already done as a first line.
• the case m == 1 unconditionally return r // n since this division is well-defined for polynomials as well.

def number_of_irreducible_polynomials(n, q=None, m=1):
    n = ZZ(n)
    if n <= 0:
        raise ValueError('n must be positive')
    if q is None:
        q = QQ['q'].gen()       # for m > 1, we produce an integer-valued polynomial in q, but it does not necessarily have integer coefficients
    R = parent(q)
    if m == 1:
        from sage.arith.misc import moebius
        r = sum((moebius(n//d) * q**d for d in n.divisors()), R.zero())
        return r // n
    elif m > 1:
        r = []
        for d in range(n):
            r.append( (q**binomial(d+m,m-1) - 1) // (q-1) * q**binomial(d+m,m) - sum(prod(binomial(r_+t-1,t) for r_,t in zip(r,p.to_exp(d))) for p in Partitions(d+1,max_part=d)) )
        return r[-1]
    else:
        raise ValueError('m must be positive')

Thanks for writing this code @maxale; I've added it to the branch together with the corresponding documentation changes.

@yyyyx4 I did something stupid in my suggestion while copy pasting and introduced whitespace issues!

Documentation preview for this PR (built with commit d952887; changes) is ready! 🎉